package ru.mersenn;

import java.util.Arrays;
import java.util.Date;

import jcuda.Pointer;
import jcuda.jcufft.JCufft;
import jcuda.jcufft.cufftHandle;
import jcuda.jcufft.cufftType;
import jcuda.runtime.JCuda;
import jcuda.runtime.cudaMemcpyKind;

import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;

/**
 * Скорость 1000 вычислений
 * q = 10^4 0.6 сек                6 сек
 * q = 10^5 6сек                   5 мин
 * q = 10^6 30 сек                 8ч
 * 
 * из статьи 57.885.161 проверяли 6 суток на 32 ядерном сервере
 * ожидается, что моя программа будет 8 * 100 = 1000 часов
 * 
 * @author akargapolov
 * @since 13.02.2013
 */
public class MersennSeq
{
	static final int ONE_DIGIT = 10;
	static final int BASE = 1 << ONE_DIGIT;

	double[][] value;
	int q;
	int realLength;
	
	int rightShift;
	int leftShift;

	public MersennSeq(int q)
	{
		realLength = q / ONE_DIGIT + 1;
		
		rightShift = q / ONE_DIGIT;
		leftShift = q % ONE_DIGIT;
		
		this.q = q;
		
		int size = 1;
		while (size <= realLength)
			size = size << 1;
		size = size << 1;
		value = new double[2][size];
		
		value[0][0] = 4;
	}

	public boolean isPrime()
	{
//		BigInteger m = new BigInteger("2").pow(q).subtract(BigInteger.ONE);
//		BigInteger l = new BigInteger("4");

	
		for (int i = 1; i <= q - 2; ++i)
		{
			if (i % 1000 == 0)
			{
				Mersenn.log.info(i);
			}
			
			pow2();
			minus2();
			modM();
			
//			l = l.pow(2).subtract(BigInteger.valueOf(2)).mod(m);
//			
//			Mersenn.log.info(l + " " + Arrays.toString(value));
		}

		return isZero();
	}

	private static long fftTime = 0;
	private void pow2()
	{
		Date start = new Date();
		
		float[] tmp = new float[value[0].length * 2];
		for (int i = 0;i<value[0].length;++i){
			tmp[2*i] = (float)value[0][i];
			tmp[2*i + 1] = 0;
		}
		float[] tmpRes = new float[tmp.length * 2];
		cufftHandle plan = new cufftHandle();
	    JCufft.cufftPlan1d(plan, tmp.length, cufftType.CUFFT_C2C, 1);
	    JCufft.cufftExecC2C(plan, tmp, tmpRes, JCufft.CUFFT_FORWARD);
	    JCufft.cufftDestroy(plan);
	     
	    double[] rrr = jcudaTransformDeviceMemory(value[0].clone());
		FastFourierTransformer.transformInPlace(value, DftNormalization.STANDARD, TransformType.FORWARD);
//		Mersenn.log.info((new Date().getTime() - start.getTime()) + "ms");
		double real, im;
		
		Mersenn.log.info(Arrays.toString(value[0]));
		Mersenn.log.info(Arrays.toString(value[1]));
//		Mersenn.log.info(Arrays.toString(tmpRes));
		Mersenn.log.info(Arrays.toString(rrr));

		
//		start = new Date();
		for (int i = 0; i < value[0].length; ++i)
		{
			real = value[0][i] * value[0][i] - value[1][i] * value[1][i];
			im = 2 * value[0][i] * value[1][i];
			
			value[0][i] = real;
			value[1][i] = im;
		}
		
		double[] mmm = new double[rrr.length];
		for (int i = 0;i<rrr.length;i+=2)
		{
			real = rrr[i] * rrr[i] - rrr[i+1]*rrr[i+1];
			im = 2 * rrr[i]*rrr[i + 1];
			
			mmm[i] = (float)real;
			mmm[i + 1] = (float)im;
		}
//		for (int i = 0;i*2<value[0].length;i++){
//			mmm[2*i] = value[0][i];
//			mmm[2*i+1]= value[1][i];
//		}
		

		plan = new cufftHandle();
		JCufft.cufftPlan1d(plan, tmpRes.length, cufftType.CUFFT_C2C, 1);
	    JCufft.cufftExecC2C(plan, tmpRes, tmpRes, JCufft.CUFFT_INVERSE);
	    JCufft.cufftDestroy(plan);
		
    
//		Mersenn.log.info((new Date().getTime() - start.getTime()) + "ms");
		FastFourierTransformer.transformInPlace(value, DftNormalization.STANDARD, TransformType.INVERSE);
		mmm = jcudaTransformDeviceMemoryInverse(mmm);
		
	    Mersenn.log.info(Arrays.toString(value[0]));
//	    Mersenn.log.info(Arrays.toString(tmpRes));
	    Mersenn.log.info(Arrays.toString(mmm));

		for (int i = 0; i < value[0].length; ++i)
		{
			value[0][i] = (long) (value[0][i] + 0.5);
			value[1][i] = 0;
		}
		
		fftTime += new Date().getTime() - start.getTime();
	}
	
	 public static double[] jcudaTransformDeviceMemory(double[] inputData){
		 
	     long timeStart = new Date().getTime();
	     double[] fftResults = new double[inputData.length + 2];
	 
	     // указатель на устройство
	     Pointer deviceDataIn = new Pointer();
	     // выделение памяти на видеокарте для входных данных
	     JCuda.cudaMalloc(deviceDataIn, inputData.length * 8);
	     // копирование входных данных в память видеокарты
	     JCuda.cudaMemcpy(deviceDataIn, Pointer.to(inputData), inputData.length * 8, 
	          cudaMemcpyKind.cudaMemcpyHostToDevice);
	      
	     Pointer deviceDataOut = new Pointer();
	     // выделение памяти на видеокарте для результатов преобразования
	     JCuda.cudaMalloc(deviceDataOut, fftResults.length * 8);
	 
	     // создание плана
	     cufftHandle plan = new cufftHandle();
	     JCufft.cufftPlan1d(plan, inputData.length, cufftType.CUFFT_D2Z, 1);
	 
	     // выполнение БПФ
	     JCufft.cufftExecD2Z(plan, deviceDataIn, deviceDataOut);
	     System.out.println("Время преобразования: " + (new Date().getTime() - timeStart)/1000.+" сек");
	      
	     // копирование результатов из памяти видеокарты в оперативную память
	     JCuda.cudaMemcpy(Pointer.to(fftResults), deviceDataOut, fftResults.length * 8, 
	             cudaMemcpyKind.cudaMemcpyDeviceToHost);
	 
	     // освобождение ресурсов
	     JCufft.cufftDestroy(plan);
	     JCuda.cudaFree(deviceDataIn);
	     JCuda.cudaFree(deviceDataOut);
	 
	  return fftResults;
	 }
	 
	 
	 public static double[] jcudaTransformDeviceMemoryInverse(double[] inputData){
		 
	     long timeStart = new Date().getTime();
	     double[] fftResults = new double[inputData.length];
	 
	     // указатель на устройство
	     Pointer deviceDataIn = new Pointer();
	     // выделение памяти на видеокарте для входных данных
	     JCuda.cudaMalloc(deviceDataIn, inputData.length * 8);
	     // копирование входных данных в память видеокарты
	     JCuda.cudaMemcpy(deviceDataIn, Pointer.to(inputData), inputData.length * 8, 
	          cudaMemcpyKind.cudaMemcpyHostToDevice);
	      
	     Pointer deviceDataOut = new Pointer();
	     // выделение памяти на видеокарте для результатов преобразования
	     JCuda.cudaMalloc(deviceDataOut, fftResults.length * 8);
	 
	     // создание плана
	     cufftHandle plan = new cufftHandle();
	     JCufft.cufftPlan1d(plan, inputData.length, cufftType.CUFFT_Z2Z, 1);
	 
	     // выполнение БПФ
	     JCufft.cufftExecZ2Z(plan, deviceDataIn, deviceDataOut, JCufft.CUFFT_INVERSE);
	     System.out.println("Время преобразования: " + (new Date().getTime() - timeStart)/1000.+" сек");
	      
	     // копирование результатов из памяти видеокарты в оперативную память
	     JCuda.cudaMemcpy(Pointer.to(fftResults), deviceDataOut, fftResults.length * 8, 
	             cudaMemcpyKind.cudaMemcpyDeviceToHost);
	 
	     // освобождение ресурсов
	     JCufft.cufftDestroy(plan);
	     JCuda.cudaFree(deviceDataIn);
	     JCuda.cudaFree(deviceDataOut);
	 
	  return fftResults;
	 }

	private void normalize()
	{
		for (int i = 0; i < value[0].length - 1; ++i)
		{
			if (value[0][i] >= BASE)
			{
				value[0][i + 1] += (long)value[0][i] / BASE;
				value[0][i] = (long)value[0][i] % BASE;
			}
			
			while (value[0][i] < 0)
			{
				if (i == (realLength - 1))
				{
					value[0][i] += (1 << leftShift);
				}
				else
				{
					value[0][i] += BASE;
					value[0][i + 1]--;
				}
			}
			
		}
	}
	
	private void minus2()
	{
		value[0][0] -= 2;
		normalize();
	}
	
	private void modM()
	{
		int i;
		boolean shifted;
		do
		{
			shifted = false;
			normalize();
			i = realLength - 1;
			if (value[0][i] >= (1 << leftShift))
			{
				shifted = true;
				value[0][i - rightShift] += ((long)value[0][i] >> leftShift);
				value[0][i] %= (1 << leftShift);
			}
			for (i = realLength; i < value[0].length; ++i)
				if (value[0][i] > 0)	//двигаем вправо на q разрядов
				{
					shifted = true;
					value[0][i - rightShift - 1] += ((long)(value[0][i] * BASE) >> leftShift);
					value[0][i] = 0;
				}
		} while (shifted);
	}

	private boolean isZero()
	{
		int i;
		for (i = 0; i < value[0].length && value[0][i] == 0.0; ++i)
			;
		
		boolean isZero = (i == value[0].length);
		
		for (i = 0; i < value[0].length; ++i)
		{
			if (i < realLength - 1 && value[0][i] != (BASE - 1))
				break;
			
			if (i == realLength - 1 && value[0][i] != ((1 << leftShift) - 1))
				break;
			
			if (i > realLength - 1 && value[0][i] != 0)
				break;
		}
		
		boolean isM = (i == value[0].length);
		
		
//		Mersenn.log.info(Arrays.toString(value));
		Mersenn.log.info("fftTime " + fftTime);
		return isZero || isM;
	}
}
